Tunable Pinning of Burst- Waves in Extended Systems with Discrete Sources 
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We study the dynamics of waves in a system of diffusively coupled discrete nonlinear sources. We 
show that the system exhibits burst waves which are periodic in a traveling-wave reference frame. 
We demonstrate that the burst waves are pinned if the diffusive coupling is below a critical value. 
When the coupling crosses the critical value the system undergoes a depinning instability via a 
saddle-node bifurcation, and the wave begins to move. We obtain the universal scaling for the mean 
wave velocity just above threshold. 
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The effect of discrete source distribution for spatially 
extended systems is a problem of interest in such dis- 
parate fields as the biophysics of the calcium release 
waves in living cells jj],§J, pinning in the dislocation mo- 
tion in crystals ||, breathers in nonlinear crystal lat- 
tices Q, Josephson junction arrays 0), and charge den- 
sity waves in one-dimensional strongly correlated electron 
systems ]|-[|. In order to elucidate the effects of dis- 
creteness, complex models of calcium release have been 
solved numerically p 4Tl|. The simple "fire-diffuse-fire" 
model constructed in ]12| displays burst (or saltatory) 
wave fronts. These fronts either propagate or they do 
not exist: they cannot undergo pinning. If a system 
with wave pinning consists of diffusively coupled discrete 
sources Q , the standard approach to its analysis has been 
to completely discretize the dynamics. This is done by re- 
placing the diffusion term with a difference scheme for the 
field at the source sites |||7|] . This simplification neglects 
the field structure between sites. The latter is crucial for 
an understanding of the system dynamics and universal 
behavior near the pinning/depinning transition. 

In this Letter we consider a discrete array of nonlin- 
ear reaction sites embedded in a continuum in which the 
reactant diffuses. We study both analytically and numer- 
ically the propagation of burst waves and their pinning, 
as well as the universal behavior of the system near the 
pinning threshold. Our model is represented by the fol- 
lowing diffusion equation with discrete nonlinear sources 
(in A^-dimcnsional space): 
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(1) 



where u is a dimensionless concentration, D is the diffu- 
sion coefficient, a is the production rate of the reactant, 
and d is the distance between neighboring sites (chan- 
nels). These sites are located at x^'s. The reaction dy- 
namics are specified by a nonlinear function f(u). To 
describe the waves of one stable phase propagating into 
another stable phase, we choose the bistable reaction dy- 
namics Jl3| . The simplest example of bistability is given 



by f(u) — —u(u — uq)(u—1) , with two stable fixed points 
u = 0, 1 , and one unstable fixed point u — uq . 

In the present Letter we study one-dimensional (ID) 
wave propagation. Then, after rescaling x — id, t — 
tor 1 , we obtain from Eq. (0) (after dropping tildes) 

u t = /3u xx + y^5(x - i)f(u) , (2) 



with the effective dimensionless diffusion coefficient /3 = 
D/ad 2 . The system dynamics are determined by the 
balance between the dissipation and the local nonlin- 
ear forcing. The latter can be expressed through its 
potential as / = —dF/du . The potential F{u) for a 
bistable system possesses two minima (stable fixed points 
of the reaction) , separated by a maximum (unstable fixed 
point of the reaction). For cubic forcing the two minima, 
u = 0, 1 , have equal energy if the potential is symmetric, 
i.e. uo = 1/2 . In this case no long-time wave propaga- 
tion is possible. The energy source for the motion is the 
difference between the depths of the potential minima, 
as determined by 7 = 1/2 — uq (for small 7 it is *~ 7). 
For 7 > , u = is a local minimum and u = 1 a global 
minimum. For 7 < , the dynamics are analogous but 
the two minima exchange their roles. We are left, there- 
fore, with two dimensionless parameters which control 
the behavior of our system: (3 and 7 . 

In Eq. (||), the effective diffusion coefficient (i is the 
only measure of how close the system is to its continuum 
limit, Ut — (3u xx + f{u) , which is realized when [3 — > 
00 . The continuous system possesses traveling- wave kink 
solutions propagating from the locally stable minimum to 
the globally stable minimum (see, e.g., |l4[). For small 7 
this solution has the form u — 1/2 [ 1 — tanh(z/2"\/2/3 ) ] + 
0(7) (in the traveling wave coordinate z = x — Ct), and 
the propagation velocity is C ~ 'Jy/P @,^|| • 

Here we investigate the dynamics for all (3, including 
the discrete limit, when |? « 1 for which the dynam- 
ics are more complicated. We have performed numer- 
ical simulations of Eq. (||). To get rid of the singular 
form of the dynamics and to ensure convergence with 
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mesh refinement, we have rewritten the equation in terms 
of a nonlocal function U(x,t) = J Q L G(x, x') u(x' , t) dx , 
where L is the system length and G is a Green-function 
of ID Laplace operator, with no-flux boundary condi- 
tion u x = at x — and Dirichlct boundary condi- 
tion u = at x = L . The Green-function has the form 
G(x, x') = (x — x 1 ) H(x — x') + x' — L with H being a 
step-function. Equation for U(x,t) was integrated im- 
plicitly, and the inverse Green-function used to restore 
u(x, i) needed to calculate the nonlinearity. 

Far from the continuum limit the system exhibits well 
developed burst waves. Figure [j] shows a space-time plot 
of the solution u{x, t) in the burst-wave regime. These 
burst waves are propagating solutions which, in the ap- 
propriate co-moving reference frame, are strictly peri- 
odic. The bursting period is uniquely determined by (3 
and 7 . 




FIG. 1. Space-time plot of u(x, t) , with black correspond- 
ing to u = 1 and white to u = . Time goes positive down- 
ward. Parameters of our model are /3 = 0.017 , 7 = 0.3 
(mo = 0.2). The system length is L — 20 , and number of 
the grid-points is 1000 . Time step dt ~ 0.007 . The period of 
the bursting is T — 35.28 , total simulation time is 3T . 

For small enough (3 the burst waves are pinned. A 
pinned solution of Eq. (0) is a static piece-wise linear 
curve connecting the stable states u = 1 and u — . The 
solution of the ID Laplace equation is linear, so the sites 
are the points where the second derivative is proportional 
to the (5-function. Thus the problem is reduced to a set 
of algebraic equations on the values of u at the sites 

& K+i + - 2ui) + f(ui) = 0, (3) 

where k, = u(xi) . 

We define the first site (going from u = 1 to u = 0) 
where the second derivative of u is positive as "the front" , 
with corresponding site number m . In the discrete limit 
(j3 <C 1) it can be shown that the values of u at the 
sites approach 1 (0) exponentially with distance from the 



front. In particular, u m+ i ~ (3 l and 1 — u m _i ~ (3 l . Let 
us consider the problem of pinning to first order in [3 . 
Then Eq. (^) at the front becomes g = u m [ — u m (uq + 
1) + uq + 2/3]- P = 0((3 2 ) . In Figure I we plot g{u m ). 
By definition of m, (u m ) xx > 0, thus it follows from 
Eq. (^) that u m < uq . Therefore, the rightmost root 
is not appropriate and the existence of solution depends 
on the value of (3 . The critical value of j3 , j3 c , at which 
stationary solutions cease to exist, corresponds to the 
situation when two remaining roots of g(u m ) , u m = u~~ 
and u m = u + merge at the graph maximum u = u max . 
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FIG. 2. A graphical representation of g(u m ) used to solve 
the stationary problem (^) to the first order in (3 . Parameter 
7 = 0.1 (uo = 0.4). Stable (unstable) front solution corre- 
sponds to u m = u~ (u + ) . 

The stability of a stationary solution of Eq. (|^), u* 
(i = 1,2,...), is determined by the eigenvalues of the 
system obtained by linearizing system (^|). This yields 
the maximum eigenvalue X m ax = 

f u (u m ) - 2(3 + 0(/3 2 ) 
which, for the roots u m = u~ and u m = u + , is negative 
and positive, respectively. Since all remaining eigenval- 
ues are negative, u = u + is a saddle and u — u~ is a 
node (see Fig. ||). If we substitute — (3 C in the expres- 
sion for Xmax we find that X ma x — [see the definition of 
g(um)]- That is, the stability boundary for the station- 
ary solution coincides with that of its existence. This 
implies that the pinning/depinning instability occurs via 
saddle-node bifurcation. 

An analysis similar to the above, but to second order 
in (3, yields the bifurcation line for the pinning/depinning 
transition in (7 , (3) space. In Figure || we compare the 
analytical bifurcation line with that obtained by direct 
numerical simulations of Eq. (^). As shown in the figure, 
the theory is in good quantitative agreement with the 
simulations, even beyond the applicability of the small (3 
perturbation theory developed above. For 7 — > 0, the 
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energy source of wave motion vanishes and consequently 
/3 C diverges. 

Let us now study the universal behavior of our system 
near the critical diffusion coefficient j3 = j3 c . First we 
note that the system dynamics are variational, i.e., 



On 



5E[u] 
Su 



(4) 



where S/Su is a variational derivative in the space of func- 
tions u{x) of a nonlinear functional E[u], given by 



E[u] = IJ(^) 2dx + Y, F ^ 



(5) 



This functional diverges for infinite systems. Therefore 
we use only its difference for different functions u{x). 
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FIG. 3. Bifurcation line f3 = /3 c (y) for pinning/depinning 
instability in the plane 7 — j3 of our model parameters. In the 
region above the line burst waves propagate and below the 
line they are pinned. 



of the unstable kink u u (x) is represented by a quadratic 
form (in the variation of the function) , having one neg- 
ative eigenvalue and a corresponding unstable direction 
n . There exist two separatrix lines of E[u] , L\ and Li , 
which start at u u and go to u s and u\ , respectively. The 
tangent directions of the separatrices at u u are collinear 
to n . The curve L12 = L\ (J L2 joins two stable minima 
u s and and contains the saddle u u . If we continue 
this curve infinitely it will contain all nodes and saddles 
oiE[u] . 

Let us introduce a normal coordinate (arclength) s on 
this curve, ds 2 = dt 2 J (du / dt) 2 dx . Figure [|(a) shows 
E(s) obtained numerically for e < . The minima and 
maxima of E(s) [see inset in Fig. 0(a)] correspond to 
stable kinks u s and unstable kinks u u , respectively. It 
can be shown, using the definitions of E and s, that the 
system dynamics on the trajectory is governed by 



ds 
It 



dE 
ds 



(6) 
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Denote a small deviation of (3 from the critical value 
as e — (3 — f3 c and consider first the case of e < . Then, 
diffusion is not strong enough to cause wave propagation, 
and there exists a linearly stable stationary kink solution 
u s (x) , corresponding to u~ in Figure 0. Let u n (x) be a 
function obtained by translating u(x) , such that u n (x) = 
u(x+n) . If u s (x) is a minimum of E[u] then u"(x) is also 
a minimum of E[u] for any integer n , implying that the 
functional E[u] has an infinite set of minima. One can 
show that E[ul(x)} - E[u s {x)} = F(u = 1) — F(u = 0) . 

Consider two adjacent minima of E[u] , u s (x) and 
ul (x) . The corresponding basins of attractions have a 
common boundary. A minimum of E[u] on the bound- 
ary is a saddle point u u (x) of the functional E[u] , which 
represents an unstable kink, corresponding to u + in Fig- 
ure ^. Due to the discrete translational symmetry of the 
system, u™(x) is also a saddle point of E[u] for any in- 
teger n. A variation of the functional E[u] in a vicinity 



FIG. 4. Dependence of the energy functional E upon the 
arclength s along the trajectory in u-space. (a) e < ; pa- 
rameters of the simulations are /3 = 0.001 , 7 = 0.3 , L = 10 , 
number of grid-points is 1000 , time step dt = 0.01 . Inset 
shows the structure of saddle-node pair, (b) e > ; parame- 
ters are the same as in Fig. |l| 

As e increases through 0, adjacent maxima and min- 
ima of E(s) merge, a saddle-node bifurcation occurs, and 
the kink begins to move. The function E(s) for this case 
is plotted in Fig. ^(b). At e = , the extrema coalesce in 
inflection points, where the first and the second deriva- 
tives are zero. The behavior of E(s) near this point is 
given by E(s) = E{0) - As 3 /3, where A > . For e > 
the expansion of E(s) acquires a linear in s part (~ e), 
so we can write E(s) = E(0) — As 3 /3 — Bes . We can 
estimate the period of the kink motion just above the 
bifurcation as 
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(7) 



where sq is a finite constant. 

For small e , T in (R) docs not depend on sq and is 
~ l/y/e. Therefore, the average kink velocity is C = 
l/T ~ near the bifurcation point. Fig. || is a plot of 
the velocity C as a function of e = (3 — (3 C obtained by 
numerical simulations of Eq. (||). The y^-scaling near 
the critical point is clearly seen. For finite e , the scaling 
breaks and, after a transient range of (3 , the expected 
continuum scaling C ~ \f]3~ appears. 

If we identify the set of points x + n as one point, a 
moving kink will be a periodic trajectory (limit cycle) on 
the reduced manifold, born by a saddle-node bifurcation. 
This type of behavior is well known for finite-dimensional 
systems JUsf . A simple exactly soluble example which 
has both the symmetries and the scaling behavior dis- 
cussed here is given by dx/dt = 1 + (1 — e) cos27ra. For 
e < the system has an infinity of fixed points and for 
e > the fixed points cease to exist and the solution 
propagates with x(t + T) = x{t) + 1 and T ~ l/V? f° r 
small e. 



were found with the velocity scaling as C ~ D/d for small 
(3 , but no pinning. The systems in ]To[ | and [ fl2| were not 
bistable which is consistent with our conjecture. In all 
models studied so far, waves fail to propagate for large 
enough site spacing. The manner in which propagation 
failure occurs differs from system to system. For example 
in the "fire-diffuse- fire" model of calcium dynamics, 
propagation failure occurs through a sequence of period- 
doubling bifurcations and crises. These facts raise an in- 
teresting open question. In what ways can waves fail to 
propagate as the source spacing is increased? To this end 
it would be useful to systematically generalize the work 
here to nonvariational systems. Such systems would be 
also more realistic as models of calcium dynamics in liv- 
ing cells. 

We are grateful to I. Aranson, S. Ponce Dawson, B. 
Ermentrout, E. Ben-Nairn, and K. Vixie for fruitful dis- 
cussions. This work was supported by the LDRD grant 
XAA7 at Los Alamos National Laboratory. 




FIG. 5. Wave propagation velocity C vs (3 — (3 C (solid 
line). Dashed lines with slope 1/2 show the scalings 
C ~ ^Je = V/3 — Pc and C ~ \f[i for discrete and continu- 
ous domain, respectively. Dashed-dotted lines illustrate the 
boundaries of continuous, discrete, and transient domains of 
(3 . The parameters are 7 = 0.1 and /3 C ~ 0.0617. 

We have shown that there exist standing kink solu- 
tions for one-component bistable reaction-diffusion sys- 
tems with discrete sources. We have also obtained the 
universal scaling for the mean velocity just above the 
pinning threshold. We conjecture that pinning is univer- 
sal in all bistable systems with discrete sources. Pinning 
was also observed in numerical simulations of a nonvari- 
ational but bistable model for the calcium release wave 
in cardiac myocytes ]LL| . These results are in contrast to 
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